import numpy as np

import math

from pyXSteam.XSteam import XSteam
from pylab import mpl

mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默认字体：解决plot不能显示中文问题
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题


def zero_pre(x):
    if x is None:
        return 0.8
    if x<0 :
        x=0
    if x>100 :
        x=100
    return x/100

def D_pre(x):
    if x is None:
        return 400
    return x/3.6

"积分定义"
class Integrator:
    def __init__(self):
        pass
    @staticmethod
    def euler(dy, y, step_size, lower_bound=float('-inf'), upper_bound=float('inf')):
        y += dy * step_size
        y = np.clip(y, lower_bound, upper_bound)
        return y

"磨煤机"
class Coal_mills:
    def __init__(self):
        self.integrator = Integrator()
    "一次风管道模型"
    def Primary_air(self,uL, uH, Wair, Tin):#一次风
        T1 = 10.3004
        T2 = 3.676511
        WLmax = 18.83
        WHmax = 55.11
        TL = 25
        TH = 300
        CL = 1.011 + TL / 25 * (1.012 - 1.011)
        CH = 1.011 + (TH - 200) / 200 * (1.028 - 1.011)
        Cin = 1.011 + (TL - 200) / 100 * (1.019 - 1.011)

        d_Wair = 1 / T1 * (-Wair + WLmax * uL + WHmax * uH)
        d_Tin = 1 / T2 * (-Tin + (CL * WLmax * uL * TL + CH * WHmax * uH * TH) / (Cin * (WLmax * uL + WHmax * uH)))
        return d_Wair, d_Tin

    "质量平衡"
    def Mass_balance(self, Wc, Wair, Tin, Mc, Mpf):
        k10 = 0.41378
        k11 = 0.070049

        Wpf = k11 * ((24.4 / 28.8) * ((273 + Tin) / 273) * ((Wair / 10)**2)) * Mpf
        d_Mc = Wc - k10 * Mc
        d_Mpf = k10 * Mc - Wpf
        return d_Mc, d_Mpf, Wpf

    def Energy_balance(self, Wc, Tin, Wair, Mc, Mpf, Tout):
        k1 = 0.0069
        k2 = 0.195494
        k3 = 0.01
        k4 = 0.001088
        k5 = 0.093383
        k6 = 0.179987
        k7 = 0.888356
        k8 = 34.206526
        k9 = 0.016556
        k12 = -0.059868
        k13 = 0.011523
        k14 = 0.262527
        k15 = 0.14986719
        Mar = 0.025
        d_Tout = (k1 * Tin + k2) * Wair + k3 * Wc - (k4 * Tout + k5) * (Wair + Wc) + k9 * (
            k6 * Mpf + k7 * Mc + k8) + k12 * Tout - k14 * k13 * Wc * Mar * Tout * (1 - math.exp(-Wair / k15))
        return d_Tout

    def simulation(self,uL, uH, Wc, Wair, Tin, Mc, Mpf, Tout, h):
        d_Wair, d_Tin = self.Primary_air(uL, uH, Wair, Tin)
        Wair = self.integrator.euler(d_Wair, Wair, h)
        Tin = self.integrator.euler(d_Tin, Tin, h)
        d_Mc, d_Mpf, Wpf = self.Mass_balance(Wc, Wair, Tin, Mc, Mpf)
        Mc = self.integrator.euler(d_Mc, Mc, h)
        Mpf = self.integrator.euler(d_Mpf, Mpf, h)
        d_Tout = self.Energy_balance(Wc, Tin, Wair, Mc, Mpf, Tout)
        Tout = self.integrator.euler(d_Tout, Tout, h)
        return Wair, Tin, Mc, Mpf, Tout, Wpf


"蒸发段"
class Evaporation_section:
    def __init__(self):
        pass
    def Evaporation_section(self,Ne, Ub, Pd, C):#蒸发段
        Ub = np.array(Ub)
        C = np.array(C)
        delta_h = 2.6665 * Pd * Pd - 83.5511 * Pd + 2055.2698 #焓差
        a = 1374.8
        b = 3397

        Dxj = (Ub@C.T*(1-0.036)-(a*Ne+b))/(delta_h*3)+180
        return Dxj

    def simulation(self, Ne, Ub, Pd, C):
        return self.Evaporation_section(Ne, Ub, Pd, C)


"省煤器"
class economizer:
    def __init__(self):
        self.integrator = Integrator()

    def economizer(self,Ne, D0, Tsmy1,Tsmj,Psm1,Tsm1,Pout):
        Mj = 5000
        cj = 0.5
        V = 80
        Vy = 100
        ky = 80

        P0 = -4.9635549e-8*Ne*Ne*Ne + 6.63318631e-5*Ne*Ne - 0.0153*Ne + 15.5287

        mid0_coeff = np.array([
            [-5.552223890598186e-5,0.188222106275285,187.4638],
            [-0.00024444444444444345,0.4033333333333324,350.0000000000002],
            [-0.0015650666666666663,2.424473333333333,-357.2819999999998],
            [-1.6666666664670333e-09,1.4033333331143116e-06,-3.699999995103936e-05],
            [-0.00018639250135555523,0.16768895880333307, -29.693576604999947]], dtype=np.float64)
        Ne_powers = np.array([[Ne ** 2], [Ne], [1]], dtype=np.float64)
        mid0 = np.dot(mid0_coeff, Ne_powers)
        T0 = mid0[0,0]
        Ty0 = mid0[1,0]
        Dy = 0.85*mid0[2,0]
        ksi1 = 20*mid0[3, 0]
        kc = 12*mid0[4, 0]

        mid1_coeff = np.array([
        [1210.86211811911,-1.59793418318532,-35.8884167229680,0.313417735621854,-0.000455489263743170,-0.867613946591206,0.00280501652204646,-0.000615644347189089],
        [-35.8884167229680,0.313417735621854,-0.867613946591206*2,2*0.00280501652204646,-0.000615644347189089,0,0,0],], dtype=np.float64)

        TP_coeff= np.array([
            [1, Tsm1, Psm1, Psm1 * Tsm1, Tsm1 ** 2, Psm1 ** 2, Psm1 ** 2 * Tsm1, Psm1 * Tsm1 ** 2]
        ], dtype=np.float64).T

        mid1 = np.dot(mid1_coeff, TP_coeff)
        rho1 = mid1[0, 0]
        drho1 = mid1[1, 0]

        H0 = 106.780612309849 + 3.54538497629407 * T0 - 16.088116133118092 * P0 + 0.063459069518316 * P0 * T0 + 0.001505531044285 * T0* T0  + 1.264454190598599 * P0 * P0 - 0.007735643634244 * P0 * P0 * T0 + 1.086607348530378e-05 * P0 * P0 * T0 * T0

        mid2_coeff = np.array([
        [2.91894501846895,4.48997695937590,-47.8864910940965,0.152989961918212,0.000109879475056696,6.18027726004829,-0.0348478482739920,4.78907792762745e-05],
        [4.48997695937590,2*0.000109879475056696,0.152989961918212,0,0,-0.0348478482739920,2*4.78907792762745e-05,0],], dtype=np.float64)

        TP_coeff1= np.array([
            [1, Tsm1, Psm1, Psm1 * Tsm1, Tsm1 ** 2, Psm1 ** 2, Psm1 ** 2 * Tsm1, Psm1 * Psm1 * Tsm1 ** 2]
        ], dtype=np.float64).T

        mid2 = np.dot(mid2_coeff, TP_coeff1)
        H1 = mid2[0, 0]
        dHT1 = mid2[1, 0]

        dHP1 = -47.8864910940965 + 0.152989961918212*Tsm1+ 2 * 6.18027726004829*Psm1 -2 * 0.0348478482739920*Psm1 * Tsm1 + 2 * 4.78907792762745e-05*Psm1 * Tsm1 * Tsm1

        mid3_coeff  = np.array([
        [1.799999999999994e-05, -0.015979999999999953, 4.03699999999999],
        [-1.9999999999988916e-07, 0.0004899999999999016, 0.9900000000000216],
        [2.2204460486140494e-20, -0.0009200000000000152, 0.8930000000000026]
        ], dtype=np.float64)

        Ty_coeff = np.array([[Tsmy1**2, Tsmy1, 1]], dtype=np.float64).T
        mid3 = np.dot(mid3_coeff, Ty_coeff)
        rhoy1 = mid3[0, 0]
        cy0 = mid3[1, 0]
        cy1 = mid3[2, 0]


        complex_num = complex(rho1 * (Psm1 - Pout) / ksi1)
        D1 = np.real(np.sqrt(complex_num))

        D1_powered = np.power(D1, 0.8)
        Dy_powered = np.power(Dy, 0.8)
        KD1 = D1_powered.real
        KDy = Dy_powered.real

        dP1 = (D0 - D1) / (V * drho1)
        dT1 = (D0 * (H0 - H1) + kc * KD1 * (Tsmj - Tsm1) - V * rho1 * dHP1 * dP1) / (V * rho1 * dHT1)
        dTj = (ky * KDy * (Tsmy1 - Tsmj) - kc * KD1 * (Tsmj - Tsm1)) / (Mj * cj)
        dTy1 = (Dy * cy0 * Ty0 - Dy * cy1 * Tsmy1 - ky * KDy * (Tsmy1 - Tsmj)) / (rhoy1 * Vy * cy1)

        return dP1, dT1, dTy1,dTj,D1

    def simulation(self, Ne, D0, Tsmy1,Tsmj,Psm1,Tsm1,Pout, h):
        dP1, dT1, dTy1,dTj,D1 = self.economizer(Ne, D0, Tsmy1,Tsmj,Psm1,Tsm1,Pout)
        Psm1 = self.integrator.euler(dP1, Psm1, h)
        Tsm1 = self.integrator.euler(dT1, Tsm1, h)
        Tsmy1 = self.integrator.euler(dTy1, Tsmy1, h)
        Tsmj = self.integrator.euler(dTj, Tsmj, h)
        return Psm1, Tsm1,Tsmy1, Tsmj,D1

"汽包"
class Drum:
    def __init__(self):
        self.a = XSteam(XSteam.UNIT_SYSTEM_MKS)
        self.integrator = Integrator()

    def Drum(self, Dsm, Pd, v, Dxj, Dgr, deta_xjs):
        R = 461.5  # 蒸汽气体常数
        vqb = 256  # 汽包
        # print("Pd received:", Pd)
        # 液体和蒸汽密度，及饱和温度的计算，使用IAPWS97库
        rho = self.a.rhoL_p(10*Pd)
        rhos = self.a.rhoV_p(10*Pd)
        T = self.a.tsat_p(10*Pd)

        # 流体和蒸汽的体积变化
        deta_V = (Dsm - Dxj) / rho
        deta_Vs = (Dxj - deta_xjs) / rhos
        deta_Pd = R * (T + 273.15) / (vqb - v) * (Dxj - Dgr)
        Ad = 45.318
        dld_dt = 1 / Ad * (deta_V + deta_Vs)

        return deta_Pd, deta_V, dld_dt, T

    def Drum_export(self, Pd, Pt):
        C0 = 26.12
        Dgr = C0 * (0.08764256902732083 * Pd * Pd - 1.8391411470804848 * Pd + 18.439682689872793) * math.sqrt(abs(Pd - Pt)) + 67.44
        return Dgr

    def simulation(self, Dsm, Pd, Pt, v, Dxj, deta_xjs, h):
        Dgr = self.Drum_export(Pd,Pt)
        deta_Pd, deta_V, dld_dt, T = self.Drum(Dsm, Pd, v, Dxj, Dgr, deta_xjs)
        Pd = self.integrator.euler(deta_Pd*0.000001, Pd, h)
        v = self.integrator.euler(deta_V, v, h)
        return Pd, T ,v,Dgr

"过热器"
class Superheater:
    def __init__(self):
        self.a = XSteam(XSteam.UNIT_SYSTEM_MKS)
        self.integrator = Integrator()
    def Superheater(self,u,Ne,P1,P2,T1,D0,T2,Tj,Ty2,Pout,P):
        Mw = 5000
        cw = 0.503
        V1 = 15
        V2 = 80
        Vt = 80
        Mj = 8000
        cj = 0.5
        Vy = 100
        ky = 1.8
        T0 = 450

        mid0_coeff = np.array([
            [2.54647189105045e-05,-0.00441307720470374,15.1093745189622],
            [-3.98703097667872e-05,0.140190285041051,119.021779144517],
            [-0.000849*0.58333338*0.58333338, 0.58333338*1.171, 780]], dtype=np.float64)
        Ne_powers = np.array([[Ne ** 2], [Ne], [1]], dtype=np.float64)
        mid0 = np.dot(mid0_coeff, Ne_powers)
        Pjw = mid0[0,0]
        Tjw = mid0[1,0]
        Ty1 = mid0[2,0]

        mid1_coeff = np.array([
        [2.85195241345586e-5,-0.0113109873432782,15.9447132979101],
        [2.5720443189266e-5, -0.0124963762893378, 15.9972804244775],
        [-0.00013200000000000234,1.2354000000000023,8.999999999999519],
        [-1.292533336e-9,9.85453335e-7,2.3517e-5],
        [-1.2679066668374196e-10,1.3287553334573863e-07,-2.218278000710685e-05],
        [-0.00016547114875555569,0.1467108342400001, -24.48438757900002]
        ], dtype=np.float64)
        mid1 = np.dot(mid1_coeff, Ne_powers)
        P0 = mid1[0,0]
        # Pout = mid1[1,0]
        Dy = mid1[2,0]
        ksi1 = mid1[3,0]
        ksi2 = mid1[4,0]
        kc = mid1[5,0]

        mid2_coeff = np.array([
        [0.00360160232767923,0.216414846272030,0.0428696651361507,-0.0163162564421020,-0.000400249182512571,0.627754112214308,-0.00111821461129329,3.89479335051162e-05],
        [0.0428696651361507,-0.0163162564421020,0.627754112214308*2,-0.00111821461129329*2,3.89479335051162e-05, 0, 0, 0],
        ], dtype=np.float64)

        TP_coeff1= np.array([
            [1, T1, P1, P1 * T1, T1 ** 2, P1 ** 2, P1 ** 2 * T1, P1 * T1 ** 2],
            [1, T2, P2, P2 * T2, T2 ** 2, P2 ** 2, P2 ** 2 * T2, P2 * T2 ** 2]
        ], dtype=np.float64).T

        mid2 = np.dot(mid2_coeff, TP_coeff1)
        rho1 = mid2[0, 0]
        rho2 = mid2[0, 1]
        drho1 = mid2[1, 0]
        drho2 = mid2[1, 1]

        mid4_coeff = np.array([
         [0.0850179093273558,11.1141641034168,0.822278858616258,-0.00343553551151210,-0.00831796694442522,10.2718063816235,-0.0385782555249321,3.62259800191877e-05],
         [11.1141641034168, -0.00831796694442522*2,-0.00343553551151210, 0, 0, -0.0385782555249321, 3.62259800191877e-05*2, 0]
         ], dtype=np.float64)
        TP_coeff2= np.array([
            [1, T1, P1, P1 * T1, T1 ** 2, P1** 2, P1 ** 2 * T1, P1 *P1* T1** 2],
            [1, T2, P2, P2 * T2, T2 ** 2, P2** 2, P2 ** 2 * T2, P2 *P2* T2** 2]
        ], dtype=np.float64).T

        mid4 = np.dot(mid4_coeff, TP_coeff2)
        H0 = 0.0850179093273558+11.1141641034168*T0+0.822278858616258*P0-0.00343553551151210*P0*T0-0.00831796694442522*T0*T0+10.2718063816235*P0*P0-0.0385782555249321*P0*P0*T0+3.62259800191877e-05*P0*P0*T0*T0
        H1 = mid4[0,0]
        H2 = mid4[0, 1]
        dHT1 = mid4[1, 0]
        dHT2 = mid4[1, 1]

        mid5_coeff = np.array([
        [0.822278858616258, -0.00343553551151210, 10.2718063816235*2,-0.0385782555249321*2, 0, 0, 0, 3.62259800191877e-05*2]
        ], dtype=np.float64)
        mid5 = np.dot(mid5_coeff, TP_coeff1)
        dHP1 = mid5[0, 0]
        dHP2 = mid5[0, 1]

        Hjw = 12.6413293982196+0.840607581456595*Pjw+3.96998415526521*Tjw-0.00575726351823391*Pjw*Pjw-0.000311550896153090*Pjw*Tjw+0.00100418117300614*Tjw*Tjw

        mid3_coeff  = np.array([
        [0.0000004, -0.00101, 0.8689],
        [-1.4444444444444416e-07, 0.00046999999999999955,0.9840000000000002]
        ], dtype=np.float64)

        Ty_coeff = np.array([[Ty1**2, Ty1, 1], [Ty2**2, Ty2, 1]], dtype=np.float64).T
        mid3 = np.dot(mid3_coeff, Ty_coeff)
        rhoy2 = mid3[0, 1]
        cy1 = mid3[1, 0]
        cy2 = mid3[1, 1]

        dPout = (D0 - P) / drho2 / Vt

        complex_num1 = complex(rho1 * (P1 - P2) / ksi1)
        complex_num2 = complex(rho2 * (P2 - Pout) / ksi2)
        D1 = np.real(np.sqrt(complex_num1))
        D2 = np.real(np.sqrt(complex_num2))

        D2_powered = np.power(D2, 0.8)
        Dy_powered = np.power(Dy, 0.8)
        KD2 = D2_powered.real
        KDy = Dy_powered.real

        dP1 = (D0 - D1) / (V1 * drho1) + u / (V1 * drho1)
        dP2 = (D1 - D2) / (V2 * drho2)
        dT1 = (D0 * (H0 - H1) - (V1 * rho1 * dHP1 - V1) * dP1) / (Mw * cw + V1 * rho1 * dHT1) + u * (Hjw - H1 + (V1 * rho1 * dHP1 - V1) / (V1 * drho1)) / (Mw * cw + V1 * rho1 * dHT1)
        dT2 = (D1 * (H1 - H2) + kc * KD2 * (Tj - T2) - V2 * rho2 * dHP2 * dP2) / (V2 * rho2 * dHT2)
        dTj = (ky * KDy * (Ty2 - Tj) - kc * KD2 * (Tj - T2)) / (Mj * cj)
        dTy2 = (Dy * cy1 * Ty1 - Dy * cy2 * Ty2 - ky * KDy * (Ty2 - Tj)) / (rhoy2 * Vy * cy2)

        return dP1, dP2, dT1, dT2, dTj, dTy2,D2,dPout

    def simulation(self, u,Ne,P1,P2,T1,D0,T2,Tj,Ty2,Pout,P, h):
        dP1, dP2, dT1, dT2, dTj, dTy2, D2, dPout = self.Superheater(u,Ne,P1,P2,T1,D0,T2,Tj,Ty2,Pout,P)
        P1 = self.integrator.euler(dP1, P1, h)
        P2 = self.integrator.euler(dP2, P2, h)
        T1 = self.integrator.euler(dT1, T1, h)
        T2 = self.integrator.euler(dT2, T2, h)
        Tj = self.integrator.euler(dTj, Tj, h)
        Ty2 = self.integrator.euler(dTy2, Ty2, h)
        Pout = self.integrator.euler(dPout, Pout, h)
        return P1, P2,T1, T2,Tj,Ty2,Pout, D2

"汽轮机做功段"
class zuogong:
    def __init__(self):
        self.integrator = Integrator()
    def zuogong(self,Ne, Ut, Pt):
        k2 = 1.1893
        D = 23.78*Pt-543.5+724.9*Ut
        dNe = (k2*D-Ne)/30
        return D, dNe

    def simulation(self, Ne, Ut, Pt, h):
        D, dNe = self.zuogong(Ne, Ut, Pt)
        Ne = self.integrator.euler(dNe, Ne, h)
        return D, Ne

"汽轮机"
class steam_turbine:

    def __init__(self,ita):
        self.a = XSteam(XSteam.UNIT_SYSTEM_MKS)
        self.integrator = Integrator()
        self.ita = ita

    "调节级"

    def regulating_level(self, Ts, Pst, Dst):
        """ 初始化"""
        # hst2=vst2=ht2=hst1=vst1=ht1=t=ht=hst=0
        """常数"""
        ptvwo = 14  # VWO工况下调节级后压力
        pstb2vwo = 16.67  # VWO工况下流经部分开启调节阀后的蒸汽压力
        Dstbvwo = 290.76  # 设计工况下流经部分开启调节阀后的蒸汽流量
        ptd = 12.3  # 设计工况下调节级后压力
        Dstd = 513.6  # 设计工况下主蒸汽流量
        fai = 1  # 喷嘴和动叶速度系数
        Kn = 1.3  # 蒸汽绝热指数
        Dstvwo = 581.5138889  # VWO工况下主蒸汽流量
        hst = self.a.h_pt(Pst*10, Ts)
        """调节级"""
        if Dst <= 290.7569444:
            Dt = Dst  # 调节级后蒸汽流量
            pt = Dst * ptd / Dstd  # 变工况下调节级后压力
            pst2 = Dst * ((pstb2vwo - ptvwo) / Dstbvwo) + pt  # 蒸汽流经部分开调节阀后的压力
            hst2 = self.a.h_pt(Pst * 10, Ts)
            vst2 = self.a.v_ph(pst2 * 10, hst2)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))
            vt2 = (pst2 / pt) ** (1 / nm) * vst2  # 通过部分开调节阀蒸汽流经喷嘴动叶后比体积
            tt2 = (1000 * pt * vt2 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过部分开调节阀蒸汽流经喷嘴动叶后温度
            ht = self.a.h_pt(pt * 10, tt2)
            t = self.a.t_ph(pt * 10, ht)
            # 调节级蒸汽温度tt
        elif 290.7569444 < Dst <= 436.1354167:
            """#部分开阀"""
            pt = Dst * ptd / Dstd  # 变工况下调节级后压力
            pst2 = (Dst - 0.5 * Dstvwo) * ((pstb2vwo - ptvwo) / Dstbvwo) + pt  # 蒸汽流经部分开调节阀后的压力
            hst2 = self.a.h_pt(Pst * 10, Ts)
            vst2 = self.a.v_ph(pst2 * 10, hst2)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))
            vt2 = (pst2 / pt) ** (1 / nm) * vst2  # 通过部分开调节阀蒸汽流经喷嘴动叶后比体积
            tt2 = (1000 * pt * vt2 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过部分开调节阀蒸汽流经喷嘴动叶后温度
            ht2 = self.a.h_pt(pt * 10, tt2)  # 通过部分开启调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值ht
            '''#全开阀'''
            hst1 = self.a.h_pt(Pst * 10, Ts)
            pst1 = 0.95 * Pst  # 蒸汽流经全开调节阀门后压力pst1
            vst1 = self.a.v_ph(pst1 * 10, hst1)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))  # 多变指数
            vt1 = (pst1 / pt) ** (1 / nm) * vst1  # 通过全开调节阀蒸汽流经喷嘴动叶后比体积
            tt1 = (1000 * pt * vt1 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过全开调节阀蒸汽流经喷嘴动叶后温度
            ht1 = self.a.h_pt(pt * 10, tt1)  # 通过全开调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值
            """%混合"""
            Dt = Dst  # 调节级后蒸汽流量
            ht = ((0.5 * Dstvwo) * ht1 + (Dst - 0.5 * Dstvwo) * ht2) / Dst  # 混合后调节级蒸汽焓值ht
            t = self.a.t_ph(pt * 10, ht)
        else:
            '''部分开阀'''
            pt = Dst * ptd / Dstd  # 变工况下调节级后压力
            pst2 = (Dst - 0.75 * Dstvwo) * ((pstb2vwo - ptvwo) / Dstbvwo) + pt  # 蒸汽流经部分开调节阀后的压力
            hst2 = self.a.h_pt(Pst * 10, Ts)  # 主蒸汽焓值
            vst2 = self.a.v_ph(pst2 * 10, hst2)  # 流经部分开阀后的蒸汽的比体积
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))
            vt2 = (pst2 / pt) ** (1 / nm) * vst2  # 通过部分开调节阀蒸汽流经喷嘴动叶后比体积
            tt2 = (1000 * pt * vt2 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过部分开调节阀蒸汽流经喷嘴动叶后温度
            ht2 = self.a.h_pt(pt * 10, tt2)  # 通过部分开启调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值ht
            # tt = IAPWS_IF97('T_ph', pt, htb - 273.15) #调节级蒸汽温度tt
            '''全开阀'''
            hst1 = self.a.h_pt(Pst * 10, Ts)  # 主蒸汽焓值
            pst1 = 0.95 * Pst  # 蒸汽流经全开调节阀门后压力pst1
            vst1 = self.a.v_ph(pst1 * 10, hst1)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))  # 多变指数
            vt1 = (pst1 / pt) ** (1 / nm) * vst1  # 通过全开调节阀蒸汽流经喷嘴动叶后比体积
            tt1 = (1000 * pt * vt1 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过全开调节阀蒸汽流经喷嘴动叶后温度
            ht1 = self.a.h_pt(pt * 10, tt1)  # 通过全开调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值
            '''混合'''
            Dt = Dst  # 调节级后蒸汽流量
            ht = ((0.75 * Dstvwo) * ht1 + (Dst - 0.75 * Dstvwo) * ht2) / Dst  # 混合后调节级蒸汽焓值ht
            t = self.a.t_ph(pt * 10, ht)
        ne = Dst*(hst-ht)
        return Dt, pt, ht, t, ne


    "压力级第一级"

    def pressure_level1(self, Dt, pt, ht, ps1):
        p1d = 6.80877688  # 设计工况下第一级抽汽压力
        xl = 0.9  # 级间效率
        Dj1d = 583.3333333  # 设计工况下第一级蒸汽流量
        p1 = p1d * Dt / Dj1d + 0.04  # 第一级抽汽压力，即第一级出口蒸汽压力
        s1 = self.a.s_ph(pt * 10, ht)  # 级组一入口蒸汽熵
        h11 = self.a.h_ps(p1 * 10, s1)  # 第一级抽汽理想焓值
        h1 = ht - xl * (ht - h11)  # 第一级实际抽汽焓值
        if ps1 > p1:
            D1 = 0
        else:
            D1 = math.sqrt(p1-ps1)*(-0.4972*p1**3+7.227*p1**2-19.44*p1+48.69)  # 第一级抽汽流量
        ne = Dt*(ht-h1)
        return D1, p1, h1, ne

    "压力级第二级"
    def pressure_level2(self, D1, p1, h1, ps2, Dt):
        p2d = 4.04  # 设计工况下第二级抽汽压力
        xl = 0.9  # 级间效率
        Dj2d = 550.548348  # 设计工况下第二级蒸汽流量
        Dj2 = Dt - D1  # 第二级蒸汽流量
        p2 = p2d * Dj2 / Dj2d  # 第二级抽汽压力，即第二级出口蒸汽压力
        s2 = self.a.s_ph(p1 * 10, h1)  # 级组二入口蒸汽熵
        h22 = self.a.h_ps(p2 * 10, s2)  # 第二级抽汽理想焓值
        h2 = h1 - xl * (h1 - h22)  # 第二级实际抽汽焓值
        if ps2 > p2:
            D2 = 0
        else:
            D2 = math.sqrt(p2-ps2)*(-2.75*p2**3+23.6*p2**2-38.75*p2+62.27)
        ne = Dj2*(h1-h2)
        return D2, p2, h2, Dj2, ne

    "再热器"

    def Reheater(self, u, Ne, P0, P1, P2, T1, D0, T2, Tj, Ty2):
        Mw = 5000
        cw = 0.503
        V1 = 15
        V2 = 80
        Mj = 8000
        cj = 0.5
        Vy = 100
        ky = 1.8
        T0 = 310

        mid0_coeff = np.array(
            [[1.87692612654277e-05, -0.00578888926835437, 8.91746543432204],
                [9.523727263669942e-07, 7.978094633458537e-04, 1.090800953098774],
                [-0.0012222222222222161, 1.8833333333333275, 459.0000000000012]], dtype = np.float64)
        Ne_powers = np.array([[Ne ** 2], [Ne], [1]], dtype = np.float64)
        mid0 = np.dot(mid0_coeff, Ne_powers)
        Pjw = mid0 [0, 0]
        Tjw = mid0 [1, 0]
        Ty1 = mid0 [2, 0]

        mid1_coeff = np.array(
            [[4.052304382533590e-06, 0.003620858815834, 0.432686520513570],
                [3.904728178104676e-06, 0.003271018799718, 0.413283907704974],
                [-8.44799999999994e-05, 0.7906559999999994, 5.760000000000159],
                [-5.899999992687796e-11, 5.8182666599523294e-08, -2.951999989046666e-06],
                [-1.851488885438405e-11, 2.1186433300572312e-08, -3.2299099949080327e-06],
                [-0.00015078811602222254, 0.13365922177666695, -22.450863535000053]], dtype = np.float64)
        mid1 = np.dot(mid1_coeff, Ne_powers)
        # P0 = mid1[0, 0]
        Pout = mid1 [1, 0]
        Dy = 5 * mid1 [2, 0]
        ksi1 = 4.231 * mid1 [3, 0]
        ksi2 = 21.15 * mid1 [4, 0]
        kc = 2 * mid1 [5, 0]

        mid2_coeff = np.array(
            [[-4.48115156278254, 0.0304999262377712, 10.2696358428991, -0.0386335598027197, -5.02653884018188e-05,
              0.581594457386798, -0.00144989003677293, 5.48743890320239e-05],
                [-0.159137508067886, 0.000531616092656617, 6.06349729532010, -0.00931510070519609,
                 -4.23506249450393e-07, 0.0784490096522479, -0.000111794576621083, 5.57847427871955e-06],
                [10.2696358428991, -0.0386335598027197, 0.581594457386798 * 2, -0.00144989003677293 * 2,
                 5.48743890320239e-05, 0, 0, 0],
                [6.06349729532010, -0.00931510070519609, 0.0784490096522479 * 2, -0.000111794576621083 * 2,
                 5.57847427871955e-06, 0, 0, 0]], dtype = np.float64)

        TP_coeff = np.array(
            [[1, T1, P1, P1 * T1, T1 ** 2, P1 ** 2, P1 ** 2 * T1, P1 ** 2 * T1 ** 2],
                [1, T2, P2, P2 * T2, T2 ** 2, P2 ** 2, P2 ** 2 * T2, P2 ** 2 * T2 ** 2]], dtype = np.float64).T

        TP_coeff1 = np.array(
            [[1, T1, P1, P1 * T1, T1 ** 2, P1 ** 2, P1 ** 2 * T1, P1 * T1 ** 2],
                [1, T2, P2, P2 * T2, T2 ** 2, P2 ** 2, P2 ** 2 * T2, P2 * T2 ** 2]], dtype = np.float64).T

        mid2 = np.dot(mid2_coeff, TP_coeff1)
        rho1 = mid2 [0, 0]
        rho2 = mid2 [1, 1]
        drho1 = mid2 [2, 0]
        drho2 = mid2 [3, 1]

        mid4_coeff = np.array(
            [[2483.42108838623, 1.89058991447985, -27.4233583398200, 0.0335060525030819, 0.000244086781910402,
              -4.83095804603545, 0.0178237649306467, -1.66247879853088e-05],
                [1.89058991447985, 0.000244086781910402 * 2, 0.0335060525030819, 0, 0, 0.0178237649306467,
                 -1.66247879853088e-05 * 2, 0]], dtype = np.float64)
        mid4 = np.dot(mid4_coeff, TP_coeff)
        H0 = 2394.59993496083 + 2.56981600432692 * T0 - 67.7424441551175 * P0 + 0.146549501898988 * P0 * T0 - 0.000987535212698787 * T0 * T0 - 16.1764009140225 * P0 * P0 + 0.0870438074160376 * P0 * P0 * T0 - 0.000123903762344842 * P0 * P0 * T0 * T0
        H1 = mid4 [0, 0]
        H2 = mid4 [0, 1]
        dHT1 = mid4 [1, 0]
        dHT2 = mid4 [1, 1]

        mid5_coeff = np.array(
            [[-27.4233583398200, 0.0335060525030819, -4.83095804603545 * 2, 0.0178237649306467 * 2, 0, 0, 0,
              -1.66247879853088e-05 * 2]], dtype = np.float64)
        mid5 = np.dot(mid5_coeff, TP_coeff1)
        dHP1 = mid5 [0, 0]
        dHP2 = mid5 [0, 1]

        Hjw = 17.1990597528719 + 0.934374962444093 * Pjw + 3.89625419335981 * Tjw - 0.00264393556446940 * Pjw * Pjw - 0.00186672781887898 * Pjw * Tjw + 0.00134403173544832 * Tjw * Tjw

        mid3_coeff = np.array(
            [[0.0000004, -0.00101, 0.8689], [-1.466666666666658e-07, 0.00047799999999999834, 0.9762666666666676]],
            dtype = np.float64)

        Ty_coeff = np.array([[Ty1 ** 2, Ty1, 1], [Ty2 ** 2, Ty2, 1]], dtype = np.float64).T
        mid3 = np.dot(mid3_coeff, Ty_coeff)
        rhoy2 = mid3 [0, 1]
        cy1 = mid3 [1, 0]
        cy2 = mid3 [1, 1]

        complex_num1 = complex(rho1 * (P1 - P2) / ksi1)
        complex_num2 = complex(rho2 * (P2 - Pout) / ksi2)
        D1 = np.real(np.sqrt(complex_num1))
        D2 = np.real(np.sqrt(complex_num2))

        D2_powered = np.power(D2, 0.8)
        Dy_powered = np.power(Dy, 0.8)
        KD2 = D2_powered.real
        KDy = Dy_powered.real

        dP1 = (D0 - D1) / (V1 * drho1) + u / (V1 * drho1)
        dP2 = (D1 - D2) / (V2 * drho2)
        dT1 = (D0 * (H0 - H1) - (V1 * rho1 * dHP1 - V1) * dP1) / (Mw * cw + V1 * rho1 * dHT1) + u * (
                Hjw - H1 + (V1 * rho1 * dHP1 - V1) / (V1 * drho1)) / (Mw * cw + V1 * rho1 * dHT1)
        dT2 = (D1 * (H1 - H2) + kc * KD2 * (Tj - T2) - V2 * rho2 * dHP2 * dP2) / (V2 * rho2 * dHT2)
        dTj = (ky * KDy * (Ty2 - Tj) - kc * KD2 * (Tj - T2)) / (Mj * cj)
        dTy2 = (Dy * cy1 * Ty1 - Dy * cy2 * Ty2 - ky * KDy * (Ty2 - Tj)) / (rhoy2 * Vy * cy2)

        return dP1, dP2, dT1, dT2, dTj, dTy2, D2

    "压力级第三级"

    def pressure_level3(self, pre, tre, Dre, ps3):
        p3d = 2.09  # 设计工况下第三级抽汽压力
        xl = 0.9  # 级间效率
        Dj3d = 520.07  # 设计工况下第三级蒸汽流量
        p3 = p3d * Dre / Dj3d - 0.01  # 第三级抽汽压力，即第三级出口蒸汽压力
        hre = self.a.h_pt(pre * 10, tre)  # 级组三入口蒸汽焓值
        s3 = self.a.s_ph(pre * 10, hre)  # 级组三入口蒸汽熵值
        h33 = self.a.h_ps(p3 * 10, s3)  # 第三级抽汽理想焓值
        h3 = hre - xl * (hre - h33)  # 第三级实际抽汽焓值
        if ps3 > p3:
            D3 = 0
        else:
            D3 = math.sqrt(p3-ps3)*(-6.841*p3**3+31.39*p3**2-17.2*p3+22.72)  # 第三级抽汽流量
        Dj3 = Dre - D3  # 第三级蒸汽流量
        ne = Dre*(hre-h3)
        return D3, p3, h3, Dj3, ne

    "压力级第四级"
    def pressure_level4(self, h3, p3, Dj3, ps4):
        p4d = 0.99  # 设计工况下第四级抽汽压力
        xl = 0.9  # 级间效率
        Dj4d = 500.07  # 设计工况下第四级蒸汽流量
        p4 = p4d * Dj3 / Dj4d - 0.04  # 第四级抽汽压力，即第四级出口蒸汽压力
        s4 = self.a.s_ph(p3 * 10, h3)  # 级组四入口蒸汽熵
        h44 = self.a.h_ps(p4 * 10, s4)  # 第四级抽汽理想焓值
        h4 = h3 - xl * (h3 - h44)  # 第四级实际抽汽焓值
        if ps4 > p4:
            D4 = 0
        else:
            D4 = math.sqrt(p4-ps4)*(-64.69 *p4**3+144.7*p4**2-47.91*p4+35.01)  # 第四级抽汽流量
        ne = Dj3*(h3-h4)
        Dj4 = Dj3 - D4  # 第四级蒸汽流量

        return D4, p4, h4, Dj4, ne

    "压力级第五级"

    def pressure_level5(self, h4, p4, Dj4, ps5):
        p5d = 0.531  # 设计工况下第五级抽汽压力
        xl = 0.9  # 级间效率
        Dj5d = 485.54  # 设计工况下第五级蒸汽流量
        p5 = p5d * Dj4 / Dj5d - 0.02  # 第五级抽汽压力，即第五级出口蒸汽压力
        s5 = self.a.s_ph(p4 * 10, h4)  # 级组五入口蒸汽熵
        h55 = self.a.h_ps(p5 * 10, s5)  # 第五级抽汽理想焓值
        h5 = h4 - xl * (h4 - h55)  # 第五级实际抽汽焓值
        if ps5 > p5:
            D5 = 0
        else:
            D5 = math.sqrt(p5-ps5)*(-365.1*p5**3+529.1*p5**2-68.05 *p5+56.83)  # 第五级抽汽流量
        # print("D5",D5)
        t5 = self.a.t_ph(p5 * 10, h5)
        ne = Dj4*(h4-h5)
        Dj5 = Dj4 - D5 # 第五级蒸汽流量

        return D5, p5, h5, Dj5, ne

    "压力级第六级"

    def pressure_level6(self, h5, p5, Dj5, ps6):
        p6d = 0.2152712 # 设计工况下第六级抽汽压力
        xl = 0.9  # 级间效率
        Dj6d = 450.07  # 设计工况下第六级蒸汽流量
        p6 = p6d * Dj5 / Dj6d - 0.0637704559442689   # 第六级抽汽压力，即第六级出口蒸汽压力
        s6 = self.a.s_ph(p5 * 10, h5)  # 级组六入口蒸汽熵
        h66 = self.a.h_ps(p6 * 10, s6)  # 第六级抽汽理想焓值
        h6 = h5 - xl * (h5 - h66)  # 第六级实际抽汽焓值
        if ps6 > p6:
            D6 = 0
        else:
            D6 = math.sqrt(p6-ps6)*(9233*p6**3-4349*p6**2+1313*p6+2.72)  # 第六级抽汽流量
        D66 = D6/30
        # print("D6",D66)
        t6 = self.a.t_ph(p6 * 10, h6)
        ne = Dj5*(h5-h6)
        Dj6 = Dj5 - D66  # 第六级蒸汽流量
        return D6, p6, h6, Dj6, ne

    "压力级第七级"

    def pressure_level7(self, h6, p6, Dj6, ps7):
        p7d = 0.0864  # 设计工况下第七级抽汽压力
        xl = 0.9  # 级间效率
        Dj7d = 420.95  # 设计工况下第七级蒸汽流量
        # Dj7d = 500
        p7 = p7d * Dj6 / Dj7d +0.0015 # 第七级抽汽压力，即第七级出口蒸汽压力
        s7 = self.a.s_ph(p6 * 10, h6)  # 级组七入口蒸汽熵
        h77 = self.a.h_ps(p7 * 10, s7)  # 第七级抽汽理想焓值
        h7 = h6 - xl * (h6 - h77)  # 第七级实际抽汽焓值
        if ps7 > p7:
            D7 = 0
        else:
            D7 = math.sqrt(p7-ps7)*(5608*p7**3-604*p7**2+2841*p7-54.62)  # 第七级抽汽流量
        t7 = self.a.t_ph(p7 * 10, h7)
        ne = Dj6*(h6-h7)
        Dj7 = Dj6 - D7  # 第七级蒸汽流量
        return D7, p7, h7, Dj7, ne

    "末级"

    def final_level(self, h7, p7, Dj7, Pc):
        xl = 0.9
        Dc = Dj7  # 末级蒸汽流量
        # sm = self.a.s_ph(p7 * 10, h7)  # 末级入口蒸汽熵
        # hm = self.a.h_ps(Pc * 10, sm)  # 末级抽汽理想焓值
        hc = self.a.hV_p(Pc * 10)
        # hc = h7 - xl * (h7 - hm) # 末级排汽焓值
        ne = Dc*(h7-hc)
        T = self.a.t_ph(Pc*10, hc) + 3
        return hc, Dc, ne, T

    # def final_level(self, h7, p7, Dj7, Pc):
    #     xl = -2.2976*(Pc/p7)**2 + 1.2718**(Pc/p7)+0.7653
    #     Dc = Dj7  # 末级蒸汽流量
    #     s7 = self.a.s_ph(p7 * 10, h7)  # 末级入口蒸汽熵
    #     hci = self.a.h_ps(Pc * 10, s7)  # 末级抽汽理想焓值
    #     # hc = self.a.hV_p(Pc * 10)
    #     hc = h7-xl*(h7-hci)
    #     ne = Dc*(h7-hc)
    #     T = self.a.t_ph(Pc*10, hc) + 3
    #     return hc, Dc, ne, T


    def simulation(self, Ts, Pst, Dst, ps1, ps2, u, Ne, P1, P2, P2_real, T1, T2, T2_real, Tj, Ty2,
                   ps3, ps4, ps5, ps6, ps7, Pc, h):
        Dt, pt, ht, t, ne0= self.regulating_level(Ts, Pst, Dst)
        D1j = Dt
        D1, p1, h1, ne1 = self.pressure_level1(Dt, pt, ht, ps1)
        D2, p2, h2, Dj2, ne2 = self.pressure_level2(D1, p1, h1, ps2, Dt)

        # dP1, dP2, dT1, dT2, dTj, dTy2, D2r = self.Reheater(u, Ne, p2, P1, P2, T1, Dj2, T2, Tj, Ty2)
        # if D2r<=200:
        #     D2r =200
        # elif D2r>=700:
        #     D2r=700
        # P1 = self.integrator.euler(dP1, P1, h)
        # P2 = self.integrator.euler(dP2, P2, h)
        # T1 = self.integrator.euler(dT1, T1, h)
        # T2 = self.integrator.euler(dT2, T2, h)
        # Tj = self.integrator.euler(dTj, Tj, h)
        # Ty2 = self.integrator.euler(dTy2, Ty2, h)

        D2r = Dj2 - D2
        # D3, p3, h3, Dj3, ne3 = self.pressure_level3(P2, T2, D2r, ps3)
        D3, p3, h3, Dj3, ne3 = self.pressure_level3(P2_real, T2_real, D2r, ps3)
        D4, p4, h4, Dj4, ne4 = self.pressure_level4(h3, p3, Dj3, ps4)
        D5, p5, h5, Dj5, ne5 = self.pressure_level5(h4, p4, Dj4, ps5)
        D6, p6, h6, Dj6, ne6 = self.pressure_level6(h5, p5, Dj5, ps6)
        D7, p7, h7, Dj7, ne7 = self.pressure_level7(h6, p6, Dj6, ps7)
        hc, Dc, ne8, Cond_T = self.final_level(h7, p7, Dj7, Pc)
        Ne = (ne0+ne1+ne2+ne3+ne4+ne5+ne6+ne7+ne8)*self.ita*0.001
        return Dt, pt, ht, D1, p1, D1j,h1, D2, p2, Dt-D1, h2, P2, T2, D2r, D3, p3, h3, Dj3, D4, p4, h4, Dj4,\
                D5, p5, h5, Dj5, D6, p6, h6, Dj6, D7, p7, h7, Dj7, hc, Dc, P1, T1, Tj, Ty2, Ne, Cond_T, ne8

"凝汽器"
class Condenser:
    def __init__(self, current_kec):
        self.a = XSteam(XSteam.UNIT_SYSTEM_MKS)
        self.integrator = Integrator
        self.current_kec = current_kec

    def Condenser(self, Pc, n, n_num, Dc, Ta):
        k = 0.04  # 翅片管束与环境空气温度的换热系数，W / (m2·℃)
        Vn = 2.5
        Ac = 1440216  # 1653379; #凝汽器总的换热面积，m2
        rhoa = 1.043  # 环境空气密度，kg / m3
        Af = 10731  # 14131; #凝汽器迎风面积，m2
        Cpa = 1.005  # 环境空气比热容，kJ / (kg·℃)
        o = 0.9  # 风机集群运行时造成每台风机流量减小的集群因子
        Qvd = 550  # 设计工况下风机空气流量，m3 / s
        nd = 110  # 设计工况下风机群转速，rpm
        Rs = 461.5  # 气体常数
        Vc = 7700 * 45  # 直接空冷凝汽器体积，m3
        p = Pc * 10
        a = XSteam(XSteam.UNIT_SYSTEM_MKS)
        Hc1 = a.hV_p(p)
        # Hc2 = a.h_pt(p, Td)
        Tc = a.tsat_p(p)
        Hlc = a.hL_p(p)
        Vf = Qvd * n_num * o / Af / nd * n
        kc = k * (np.power(complex(Vf / Vn), 0.6)).real
        NTU = kc * Ac / rhoa / Vf / Af / Cpa
        Qc = rhoa * Vf * Af * Cpa * (1 - math.exp(-NTU)) * (Tc - Ta)
        Dec = Qc / (Hc1 - Hlc)* self.current_kec
        deta_Pc = Rs * (Tc + 273.15) * (Dc - Dec) / Vc
        # Tning =Dc * (Hc-Hlc) / (Af * Vf * rhoa * Cpa) / (1 - math.exp(-NTU)) + Ta - 8
        # Pc = 9.81 * 0.001 * 0.001 * ((Tning+100)/57.66) ** 7.46
        # Pc = 5.3699 *0.001 * ((Tning + 49.6043) / 86.0686) ** 5.3699 - 0.0075 + 0.0008
        # print("Pc",Pc)
        return deta_Pc, Dec, Tc, Hc1
        # return Pc, Dec, Tc, Vf, Hc1, Hc2
        # return Pc, Dec, Tc, Hc1

    def simulation(self, Pc, n, n_num, Dc, Ta, h):
        deta_Pc, Dec, Tc, Hc1 = self.Condenser(Pc, n, n_num, Dc, Ta)
        # Pc, Dec, Tc,  Hc1,  = self.Condenser(Hc, Dc, Pc, n, Ta)
        Pc = self.integrator.euler(deta_Pc*0.000001, Pc, h)
        return Pc, Dec, Tc, Hc1

"高低加除氧器"
class Heat_recovery_system:
    def __init__(self):
        self.a = XSteam(XSteam.UNIT_SYSTEM_MKS)
        self.integrator = Integrator()

    def Deaerator(self, Dfw,Dw5,Dd3,hd3,D4,h4,Dst,Ps4,Tw5,drhol4_dt):
        P = Ps4 * 10  # 转换单位Ps6:MPa  P：bar
        hl4 = self.a.hL_p(P)  # 除氧器压力 ps4下对应的饱和水焓，kJ/kg
        hv4 = self.a.hV_p(P)
        cp = self.a.CpL_p(P)  ##给水的定压比热容，kJ/(kg·℃)
        Vv4 = 110  # 除氧器内饱和蒸汽的体积，m3；
        Vl4 = 245  # 除氧器内饱和水体积，m3
        Ts4 = self.a.tsat_p(P)  # 第5号回热器压力下的饱和温度，℃
        detaTs4 = 3.11278e-5 * Dst * Dst - 0.02412 * Dst + 8.196209
        qcd4 = Dw5 * (hl4 - cp * Tw5) + Dd3 * (hl4 - hd3) + D4 * (hv4 - h4)
        De4 = qcd4 / (hv4 - hl4)
        deta_Ps4 = 461.5 * (Ts4 + 273.15) / (Vv4 + Vl4) * (D4 - De4)
        Tw4 = Ts4 + 0.7 * detaTs4
        # 除氧器水位
        Ld=31.656#除氧器长度
        r=1.9#除氧器半径
        hw5 = self.a.h_pt(P, Tw5)
        hs4 = self.a.hL_p(P)
        rhol4 = self.a.rhoL_p(P)#除氧器内饱和水的密度，kg/m3
        # Vcy = math.pi * Ls * Ls * (3 * r - Ls) / 3 + (0.5 * (0.5 * math.pi + math.asin((Ls - r) / r)) * r * r + (Ls - r) * np.sqrt(r * r - (r - Ls) ** 2)) * Ld
        # dLs_dt = (D4 + Dd3 + Dw5 - Dfw - Vcy * drhol4_dt) / rhol4 * 1 / (math.pi / 3 * (6 * Ls * r - 3 * Ls * Ls) + Ld * (0.5 * r / np.sqrt(1 - (Ls - r) ** 2 / r ** 2) + np.sqrt(r * r - (r - Ls) ** 2) - (Ls - r) ** 2 / np.sqrt(r * r - (r - Ls)**2)))
        return deta_Ps4, Tw4, hw5, hs4, rhol4
    def pressure_heater(self, Dfw, V3, V31, V32, V33, m33, m32, m31, ml3, kcd3, kdr3, drhodt33, drhodt32, drhodt31, Tw4, Tn3, deta_Tw4, Tl3, Tw3, D3, h3, hd3,  Dd2, hd2, Ps3):
        P = Ps3 * 10  # 转换单位Ps: MPa P：bar
        td3 = self.a.t_ph(P, hd3)
        ts3 = self.a.tsat_p(P)  # 第3号回热器压力下的饱和温度，℃
        hv3 = self.a.hV_p(P)
        hl3 = self.a.hL_p(P)
        cp = self.a.CpL_p(P)
        qdr3 = kdr3 * (np.power(complex(Dfw), 0.8)).real * 0.5 * (ts3 + td3 - Tn3 - Tw4)
        deta_Tn3 = (2 * (Dfw * (Tw4 - Tn3) + qdr3 / cp) - (Tw4 + Tn3) * drhodt33 * V33 * deta_Tw4) / m33 - deta_Tw4
        qcd3 = kcd3 * (np.power(complex(Dfw), 0.8)).real * (ts3 - (Tn3 + Tl3) / 2)
        deta_Tl3 = (2 * (Dfw * (Tn3 - Tl3) + qcd3 / cp) - (Tn3 + Tl3) * drhodt32 * V32 * deta_Tn3) / m32 - deta_Tn3
        qgr3 = D3 * (h3 - hl3)
        deta_Tw3 = (2 * (Dfw * (Tl3 - Tw3) + qgr3 /cp) - (Tl3 + Tw3) * drhodt31 * V31 * deta_Tl3) / m31 - deta_Tl3
        De3 = qcd3 / (hv3 - hl3)
        Dd3 = De3 + Dd2
        deta_ps3 = 461.5 * (ts3 + 273.15) / V3 * (D3 - De3)
        deta_hd3 = (De3 * hl3 - Dd3 * hd3 + Dd2 * hd2 - qdr3) / ml3

        return deta_Tn3, deta_Tl3, deta_Tw3, deta_ps3, Dd3, deta_hd3

    def simulation(self, Twc, Tn7, Tl7, Ps7, Tw7, D7, h7, hd7, Dd6, hd6,
                   Tn6, Tl6, Ps6, Tw6, D6, h6, Dd5, hd5,
                   Tn5, Tl5, Ps5, Tw5, D5, h5,
                   Ps4, Dd3, hd3, D4, h4, D0,
                   Tn3, Tl3, Ps3, Tw3, D3, h3, Dd2, hd2,
                   Tn2, Tl2, Ps2, Tw2, D2, h2, Dd1, hd1,
                   Tn1, Tl1, Ps1, Tw1, D1, h1,
                   Dwc_real, Dfw_real, Tw4,Twc_m,drhol4_dt,rhol4_m, h
                   ):
        ml7 = 400
        m73 = 1247
        m72 = 1247
        m71 = 1247
        V7 = 80
        V73 = 1.3  # 第7号回热加热器疏水冷却段给水体积，m3
        V72 = 1.3  # 第7号回热加热器饱和蒸汽凝结段给水体积，m3
        V71 = 1.3
        drhodt73 = -0.0053431932 * Twc - 0.1954681
        drhodt72 = -0.0053431932 * Tn7 - 0.1954681
        drhodt71 = -0.0053431932 * Tl7 - 0.1954681
        # kcd7 = 43.0445 * Ps7 + 2.7598
        kcd7 = 6.6
        kdr7 = 78.0536 * Ps7 + 2.9878
        deta_Twc = (Twc-Twc_m)*h
        Twc_m = Twc
        deta_Tn7, deta_Tl7, deta_Tw7, deta_Ps7, Dd7, deta_hd7 = self.pressure_heater(
                Dwc_real * 0.5, V7, V71, V72, V73, m73, m72, m71, ml7, kcd7, kdr7, drhodt73, drhodt72, drhodt71, Twc,
                Tn7, deta_Twc, Tl7, Tw7, D7, h7, hd7, Dd6 * 0.5, hd6, Ps7)
        Ps7 = self.integrator.euler(deta_Ps7, Ps7 * 1000000, h, lower_bound = 0, upper_bound = 10000000) / 1000000
        hd7 = self.integrator.euler(deta_hd7, hd7, h, lower_bound = 0, upper_bound = 5000)
        Tw7 = self.integrator.euler(deta_Tw7, Tw7, h, lower_bound = 0, upper_bound = 1000)
        Tl7 = self.integrator.euler(deta_Tl7, Tl7, h, lower_bound = 0, upper_bound = 1000)
        Tn7 = self.integrator.euler(deta_Tn7, Tn7, h, lower_bound = 0, upper_bound = 1000)


        ml6 = 956.3157
        m61 = 2457.6
        m62 = 2457.6
        m63 = 2457.6
        V6 = 43.83
        V61 = 2.6
        V62 = 2.6
        V63 = 2.6  # 第6号回热加热器疏水冷却段给水体积，m3
        drhodt63 = -0.004546250 * Tw7 - 0.26483375
        drhodt62 = -0.004546250 * Tn6 - 0.26483375
        drhodt61 = -0.004546250 * Tl6 - 0.26483375
        # kcd6 = 9.6606 * Ps6 + 4.2743 + 17
        kcd6 = -678.0962 * Ps6 + 112.9011
        kdr6 = 14.2294 * Ps6 + 5.5014
        deta_Tn6, deta_Tl6, deta_Tw6, deta_Ps6, Dd6, deta_hd6 = self.pressure_heater(
                Dfw_real, V6, V61, V62, V63, m63, m62, m61, ml6, kcd6, kdr6, drhodt63, drhodt62, drhodt61, Tw7, Tn6, deta_Tw7,
                Tl6, Tw6, D6, h6, hd6, Dd5, hd5, Ps6)
        Ps6 = 0.000001 * self.integrator.euler(deta_Ps6, Ps6 * 1000000, h, lower_bound = 0, upper_bound = 10000000)
        hd6 = self.integrator.euler(deta_hd6, hd6, h, lower_bound = 0, upper_bound = 5000)
        Tw6 = self.integrator.euler(deta_Tw6, Tw6, h, lower_bound = 0, upper_bound = 1000)
        Tn6 = self.integrator.euler(deta_Tn6, Tn6, h, lower_bound = 0, upper_bound = 1000)
        Tl6 = self.integrator.euler(deta_Tl6, Tl6, h, lower_bound = 0, upper_bound = 1000)


        ml5 = 1041.052632
        m51 = 2406
        m52 = 2406
        m53 = 2406
        kcd5 = 3.8604 * Ps5 + 4.3508 + 5.35  # 第5号回热加热器饱和蒸汽凝结段的传热系数(含面积)，W /℃
        kdr5 = 5.742 * Ps5 + 5.6127  # 第5号回热加热器疏水冷却段的传热系数(含面积)，W /℃
        V5 = 43.83
        V53 = 2.6  # 第5号回热加热器疏水冷却段给水体积，m3
        V52 = 2.6  # 第5号回热加热器饱和蒸汽凝结段给水体积，m3
        V51 = 2.6
        drhodt53 = -0.004423 * Tw6 - 0.276974
        drhodt52 = -0.004423 * Tn5 - 0.276974
        drhodt51 = -0.004423 * Tl5 - 0.276974

        deta_Tn5, deta_Tl5, deta_Tw5, deta_Ps5, Dd5, deta_hd5 = self.pressure_heater(
                Dfw_real, V5, V51, V52, V53, m53, m52, m51, ml5, kcd5, kdr5, drhodt53, drhodt52, drhodt51, Tw6, Tn5, deta_Tw6,
                Tl5, Tw5, D5, h5, hd5, 0, 0, Ps5)
        Ps5 = 0.000001 * self.integrator.euler(deta_Ps5, Ps5 * 1000000, h, lower_bound = 0, upper_bound = 10000000)
        hd5 = self.integrator.euler(deta_hd5, hd5, h, lower_bound = 0, upper_bound = 5000)
        Tw5 = self.integrator.euler(deta_Tw5, Tw5, h, lower_bound = 0, upper_bound = 1000)
        Tn5 = self.integrator.euler(deta_Tn5, Tn5, h, lower_bound = 0, upper_bound = 1000)
        Tl5 = self.integrator.euler(deta_Tl5, Tl5, h, lower_bound = 0, upper_bound = 1000)

        Tw4_m = Tw4

        deta_Ps4, Tw4, hw5, hs4, rhol4 = self.Deaerator(
                Dfw_real, Dwc_real, Dd3, hd3, D4, h4, D0, Ps4, Tw5, drhol4_dt)
        Ps4 = 0.000001 * self.integrator.euler(deta_Ps4, Ps4 * 1000000, h, lower_bound = 0, upper_bound = 10000000)
        # Ls = self.integrator.euler(dLs_dt, Ls, h)

        drhol4_dt = (rhol4 - rhol4_m)*h
        rhol4_m = rhol4
        deta_Tw4 = (Tw4 - Tw4_m)*h

        V3 = 38
        V31 = 0.61578947
        V32 = 6.568421053
        V33 = 0.61578947
        m33 = 334.5918679
        m32 = 5176.010333
        m31 = 321.4768759
        ml3 = 2692.1676
        kcd3 = 1.2717 * Ps3 + 5.1268 + 3.75
        kdr3 = 2.0352 * Ps3 + 4.8479
        drhodt33 = -0.009026 * Tw3 + 0.8595
        drhodt32 = -0.009026 * Tn3 + 0.8595
        drhodt31 = -0.009026 * Tl3 + 0.8595
        deta_Tn3, deta_Tl3, deta_Tw3, deta_ps3, Dd3, deta_hd3 = self.pressure_heater(
                Dfw_real, V3, V31, V32, V33, m33, m32, m31, ml3, kcd3, kdr3, drhodt33, drhodt32, drhodt31, Tw4+3.6, Tn3,
                deta_Tw4, Tl3, Tw3, D3, h3, hd3, Dd2, hd2, Ps3)
        Ps3 = 0.000001 * self.integrator.euler(deta_ps3, Ps3 * 1000000, h, lower_bound = 0, upper_bound = 10000000)
        hd3 = self.integrator.euler(deta_hd3, hd3, h, lower_bound = 0, upper_bound = 5000)
        Tw3 = self.integrator.euler(deta_Tw3, Tw3, h, lower_bound = 0, upper_bound = 1000)
        Tn3 = self.integrator.euler(deta_Tn3, Tn3, h, lower_bound = 0, upper_bound = 1000)
        Tl3 = self.integrator.euler(deta_Tl3, Tl3, h, lower_bound = 0, upper_bound = 1000)

        V2 = 37
        V21 = 0.61578947
        V22 = 6.568421053
        V23 = 0.61578947
        m23 = 355.2263415
        m22 = 5331.756671
        m21 = 337.4301625
        ml2 = 2471.111503
        kcd2 = 0.7816 * Ps2 + 5.1264 + 1.35
        kdr2 = 1.1847 * Ps2 + 6.1126
        drhodt23 = -0.0062348 * Tw2 + 0.16187
        drhodt22 = -0.0062348 * Tn2 + 0.16187
        drhodt21 = -0.0062348 * Tl2 + 0.16187
        deta_Tn2, deta_Tl2, deta_Tw2, deta_ps2, Dd2, deta_hd2 = self.pressure_heater(
                Dfw_real, V2, V21, V22, V23, m23, m22, m21, ml2, kcd2, kdr2, drhodt23, drhodt22, drhodt21, Tw3, Tn2,
                deta_Tw3, Tl2, Tw2, D2, h2, hd2, Dd1, hd1, Ps2)
        Ps2 = 0.000001 * self.integrator.euler(deta_ps2, Ps2 * 1000000, h, lower_bound = 0, upper_bound = 10000000)
        Tn2 = self.integrator.euler(deta_Tn2, Tn2, h, lower_bound = 0, upper_bound = 1000)
        Tl2 = self.integrator.euler(deta_Tl2, Tl2, h, lower_bound = 0, upper_bound = 1000)
        Tw2 = self.integrator.euler(deta_Tw2, Tw2, h, lower_bound = 0, upper_bound = 1000)
        hd2 = self.integrator.euler(deta_hd2, hd2, h, lower_bound = 0, upper_bound = 5000)

        V1 = 34
        V11 = 0.61578947
        V12 = 6.568421053
        V13 = 0.61578947
        m13 = 368.1943241
        m12 = 5822.451995
        m11 = 356.8684455
        ml1 = 2146.214863
        kcd1 = 0.6091 * Ps1 + 5.8856 + 1.75
        kdr1 = 0.9191 * Ps1 + 6.2908
        drhodt13 = -0.009026 * Tw1 + 0.8595
        drhodt12 = -0.009026 * Tn1 + 0.8595
        drhodt11 = -0.009026 * Tl1 + 0.8595
        deta_Tn1, deta_Tl1, deta_Tw1, deta_ps1, Dd1, deta_hd1 = self.pressure_heater(
                Dfw_real, V1, V11, V12, V13, m13, m12, m11, ml1, kcd1, kdr1, drhodt13, drhodt12, drhodt11, Tw2, Tn1,
                deta_Tw2, Tl1, Tw1, D1, h1, hd1, 0, 0, Ps1)
        Ps1 = 0.000001 * self.integrator.euler(deta_ps1, Ps1 * 1000000, h, lower_bound = 0, upper_bound = 10000000)
        Tn1 = self.integrator.euler(deta_Tn1, Tn1, h, lower_bound = 0, upper_bound = 1000)
        Tl1 = self.integrator.euler(deta_Tl1, Tl1, h, lower_bound = 0, upper_bound = 1000)
        Tw1 = self.integrator.euler(deta_Tw1, Tw1, h, lower_bound = 0, upper_bound = 1000)
        hd1 = self.integrator.euler(deta_hd1, hd1, h, lower_bound = 0, upper_bound = 5000)

        return Ps7, Dd7*2, hd7, Tn7, Tl7, Tw7, Ps6, Dd6, hd6, Tn6, Tl6, Tw6, Ps5, Dd5, hd5, Tn5, Tl5, Tw5,\
                Ps4, Tw4, Ps3, Dd3, hd3, Tn3, Tl3, Tw3, Ps2, Dd2, hd2, Tn2, Tl2, Tw2, Ps1, Dd1, hd1, Tn1, Tl1, Tw1,\
                rhol4,Twc_m,Tw4_m,rhol4_m,drhol4_dt

